# load and view data
load("bc_bear_occurrences.Rda")
# str(occ_data)
load("BC_Covariates.Rda")
summary(DATA)
##            Length Class           Mode
## Window      1     SpatialPolygons S4  
## Elevation  10     im              list
## Forest     10     im              list
## HFI        10     im              list
## Dist_Water 10     im              list

Visualize bear occurrences in BC

# extract location columns
bears_loc <- occ_data[, c("decimalLongitude", "decimalLatitude", "month", "year")]
bears_loc_filtered <- subset(bears_loc, year %in% c(2020, 2021, 2022, 2023, 2024))

# create sf object
bears_sf <- st_as_sf(
  bears_loc_filtered,
  coords = c("decimalLongitude", "decimalLatitude"),
  crs = 4326  # WGS84 (longitude/latitude)
)

# BC Albers projection
bears_sf_proj <- st_transform(bears_sf, crs = 3005)

# extract BC window
window_sf <- st_as_sf(DATA$Window) # convert SpatialPolygons to Simple Features (sf)
window_proj <- st_transform(window_sf, crs = 3005) # ensure same CRS
window <- as.owin(window_proj) # convert to owin using sf object

# To avoid plotting points outside window (illegal points)
# Intersect points with window
bears_sf_proj <- subset(bears_sf_proj, st_within(bears_sf_proj, window_proj, sparse = FALSE))

# extract x, y coordinates
coords <- st_coordinates(bears_sf_proj)

# create ppp object
bears_ppp <- ppp(
  x = coords[,1],
  y = coords[,2],
  window = window
)
## Warning: data contain duplicated points
plot(bears_ppp, pch = 21, main = "Black bear occurrences in BC, 2020-2024")

Compare seasonal distributions

# Define the seasons
seasons <- list(
  winter = c(12, 1, 2),
  spring = c(3, 4, 5),
  summer = c(6, 7, 8),
  autumn = c(9, 10, 11)
)

# Create an empty list to store the ppp objects for each season
ppp_list <- list()

# Create an empty list to store the filtered data for each season
bears_sf_list <- list()

# Loop over seasons
for (i in 1:length(seasons)) {
  
  # Filter for each season
  season_name <- names(seasons)[i]
  season_months <- seasons[[i]]
  
  # Filter bears_loc_filtered by the season's months
  bears_loc_season <- bears_loc_filtered[bears_loc_filtered$month %in% season_months, ]
  
  # Print the number of observations for this season
  cat(season_name, ": ", nrow(bears_loc_season), " observations\n", sep="")
  
  # Create sf object for the filtered season
  bears_sf_season <- st_as_sf(
    bears_loc_season,
    coords = c("decimalLongitude", "decimalLatitude"),
    crs = 4326  # WGS84 (longitude/latitude)
  )
  
  # Transform to BC Albers projection
  bears_sf_season_proj <- st_transform(bears_sf_season, crs = 3005)
  
  # Store sf object in the list
  bears_sf_list[[season_name]] <- bears_sf_season_proj
  
  # Extract x, y coordinates
  coords <- st_coordinates(bears_sf_season_proj)
  
  # Create ppp object for each season
  suppressWarnings(
    bears_ppp_season <- ppp(
      x = coords[, 1],
      y = coords[, 2],
      window = window
    )
  )
  
  # Store ppp object in the list
  ppp_list[[season_name]] <- bears_ppp_season
}
## winter: 76 observations
## spring: 808 observations
## summer: 1862 observations
## autumn: 841 observations
# Plot the four maps in a 2x2 grid
par(mfrow=c(2,2), mar=c(1,1,1,1))  # Set up a 2x2 plotting window

# Loop over ppp objects and plot them
for (i in 1:length(ppp_list)) {
  season_name <- names(ppp_list)[i]
  suppressWarnings(
    plot(ppp_list[[season_name]], main = season_name, pch = 21)
  )
}

# Visualize covariates
par(mfrow=c(2,2), mar=c(1,1,1,1))  # Set up a 2x2 plotting window
plot(DATA$Elevation)
plot(DATA$Forest)
plot(DATA$HFI)
plot(DATA$Dist_Water)

# Isolate covariates
elev <- DATA$Elevation
cover <- DATA$Forest
dist_water <- DATA$Dist_Water
hfi <- DATA$HFI

Estimate intensity of bear occurrences based on forest cover

# intensity as a function of forest cover overall
rho <- rhohat(bears_ppp, cover) 
plot(rho, 
     xlim=c(0, max(cover)), 
     main="Estimated Rho vs. Elevation")

# intensity as a function of forest cover broken up by seasons
par(mfrow=c(2,2), mar=c(2,2,2,2))

for (i in 1:length(ppp_list)) {
  season_name <- names(ppp_list)[i]
  rho <- rhohat(ppp_list[[season_name]], cover)
  plot(rho, 
       xlim=c(0, max(cover)), 
       main = paste(season_name))
}

# convert im to raster
im2raster <- function(im_obj) {
  m <- as.matrix(im_obj) # image to matrix
  r <- raster(m) # create raster object
  
  # range and extension
  ext <- c(im_obj$xrange[1] - im_obj$xstep/2,
           im_obj$xrange[2] + im_obj$xstep/2,
           im_obj$yrange[1] - im_obj$ystep/2,
           im_obj$yrange[2] + im_obj$ystep/2)
  extent(r) <- ext
  return(r)
}

# extract forest cover values at occurrence locations
cover_raster <- im2raster(cover)
bears_sf_proj$cover_value <- extract(cover_raster, bears_sf_proj)

KDE to estimate the concentration of bear occurrences based on forest cover values

# KDE for Bear Locations
kde_b <- density(bears_sf_proj$cover_value, na.rm = TRUE, bw = "SJ-dpi")
summary(bears_sf_proj$cover_value)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   25.44   52.25   53.14   80.95  100.00    1444
plot(kde_b,
     main = "Kernel Density Estimate of Forest Cover at Bear Locations",
     xlab = "Elevation (meters)",
     ylab = "Density",
     col = "#2C6B2F", 
     lwd = 2)

These results suggest that bears are found in a mix of forested and non-forested areas, but on average, they’re in areas with moderate to high forest cover.

Model Forest Cover

# NA with mean forest cover
mean_cover <- mean(as.vector(as.matrix(cover)), na.rm = TRUE)
cover_clean <- eval.im(ifelse(is.na(cover), mean_cover, cover))

# polynomial model
cover_model <- ppm(bears_ppp ~ cover + I(cover^2), covariates = list(cover = cover_clean))
print(cover_model)
## Nonstationary Poisson process
## Fitted to point pattern dataset 'bears_ppp'
## 
## Log intensity:  ~cover + I(cover^2)
## 
## Fitted trend coefficients:
##   (Intercept)         cover    I(cover^2) 
## -2.034720e+01  6.209799e-02 -6.168772e-04 
## 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -2.034720e+01 4.834276e-02 -2.044195e+01 -2.025245e+01   ***
## cover        6.209799e-02 2.093576e-03  5.799466e-02  6.620132e-02   ***
## I(cover^2)  -6.168772e-04 2.010334e-05 -6.562791e-04 -5.774754e-04   ***
##                   Zval
## (Intercept) -420.89447
## cover         29.66120
## I(cover^2)   -30.68531
# predicts and plots intensity
pred_intensity_cover <- predict(cover_model)
plot(pred_intensity_cover, 
     main = "Predicted Intensity Based on Forest Cover")
points(bears_ppp, 
       pch = 20, 
       col = "red")

The positive coefficient on the linear term indicates that bear occurrence intensity initially increases with forest cover. However, the negative coefficient on the quadratic term means that this relationship has a peak, after which further increases in forest cover are associated with a decrease in expected intensity. This suggests that bears are most likely to be found in areas with moderate forest cover, and less likely to occur in areas with either very low or very high forest density.

This pattern may reflect ecological preferences, where bears favor mixed habitats over open or densely forested areas. It’s also possible that human detection rates are lower in areas of dense vegetation, where visibility is reduced.

Compute and plot residuals

cover_residuals <- residuals(cover_model, type = "pearson")

# remove any non-finite residuals
cover_residuals$v[!is.finite(cover_residuals$v)] <- NA

# plot forest cover residuals
plot(cover_residuals, 
     main = "Residuals - Forest Cover Model", 
     na.col = "transparent")

Partial Residuals

# Computes partial residuals for forest cover
par_res_cover_all <- parres(cover_model, "cover")

# Plots partial residuals
plot(par_res_cover_all,
     legend = FALSE,
     lwd = 2,
     main = "Partial Residuals Forest Cover - All Seasons",
     xlab = "Elevation (m)",
     ylab = "Partial Residuals")

Seasonal Models

models_cover_seasonal <- list()

# model for each season (individually)
for (season in names(ppp_list)) {
  models_cover_seasonal[[season]] <- ppm(ppp_list[[season]], ~ cover + I(cover^2), covariates = list(cover = cover_clean))
  # cat("Model for", season, ":\n")
  # print(models_seasonal[[season]])
  plot(predict(models_cover_seasonal[[season]]), main = paste("Intensity ~ Forest Cover -", season))
}
## Warning: glm.fit: algorithm did not converge

Winter: Winter: Intercept: -23.78; Forest Coefficient: 0.0582; Quadratic Coefficient: -0.000672.

Bear occurrences increase with forest cover up to a point, then decline. The smaller linear coefficient suggests a slower rise in intensity compared to other seasons. The model did not converge, so caution is needed in interpreting the results.

Spring: Intercept: -21.96; Forest Coefficient: 0.0657; Quadratic Coefficient: -0.000632.

The relationship is similar to winter but slightly stronger, meaning a more noticeable increase and peak in moderate forest areas.

Summer: Intercept: -21.05; Forest Coefficient: 0.0621; Quadratic Coefficient: -0.000602.

Summer shows a comparable trend with a slightly weaker curvature, suggesting bears may tolerate higher forest density before intensity starts to drop.

Autumn: Intercept: -21.70; Forest Coefficient: 0.0678; Quadratic Coefficient: -0.000732.

Autumn has the sharpest curve (greatest linear and quadratic coefficient magnitudes), indicating a more distinct peak and a stronger decline in high-density forest areas.

Overall, all seasons display a consistent non-linear relationship where bear occurrence intensity increases with forest cover up to a moderate level, then declines. While the shape of the curve remains similar, the strength and location of the peak vary slightly by season, indicating subtle shifts in habitat preference throughout the year.

Partial Residuals per Season

# 2x2 plot layout
par(mfrow = c(2, 2))

# Generates and plots partial residuals per season
for (season in names(models_cover_seasonal)) {
  model <- models_cover_seasonal[[season]]
  parres_season <- parres(model, "cover")
  plot(parres_season,
       legend = FALSE,
       lwd = 2,
       main = paste("Partial Residuals - Forest Cover", season),
       xlab = "Elevation (m)",
       ylab = "Partial Residuals")
}

KDE for each season

# KDE comparison
dens_list <- lapply(ppp_list, function(p) density(p, sigma = 10000))
for (s in names(dens_list)) {
  plot(dens_list[[s]], main = paste("KDE Surface -", s))
}

Performing Quadrat test for each season

# quadrat tests
for (s in names(ppp_list)) {
  qt <- quadrat.test(ppp_list[[s]], nx = 5, ny = 5)
  print(qt)
  plot(qt, main = paste("Quadrat test -", s))
}
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 561.39, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 2202.7, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 3767, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 2720.5, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)

All p-values are small (< 2.2e-16) indicating significant deviation from CSR in each season.